Genomics, Proteomics & Bioinformatics

Tuberculosis remains a major global public health threat. While whole-genome sequencing has transformed our understanding of the causative agent, Mycobacterium tuberculosis (MTB), existing genomic databases are highly fragmented and often underrepresent structural variations (SVs). Furthermore, critical population-genetic statistics are rarely integrated with phylogenetic and geographic context, forcing researchers to reconcile separate datasets manually. To address this gap, we developed TBpop (https://tbpop.chinacdc.cn), an open-access, integrated population genomics portal. TBpop is built from 420 clinical MTB isolates selected from the first national drug resistance baseline survey in China. The portal integrates isolate metadata, pangenome categories, SNPs, SVs, IS6110 insertion sites, strain phylogeny, and gene-level statistics, and provides three interactive explorer modules: the Population Explorer, the Statistics Explorer, and the Variation Explorer. Additionally, a User Analysis module allows researchers to run population genetic workflows on their own alignments. TBpop provides an integrated platform for exploring genome plasticity, signatures of positive selection, and conservation patterns of functionally important genes in MTB.

Hepatocellular Carcinoma (HCC) arising from Metabolic Dysfunction-Associated Steatotic Liver Disease (MASLD) is an increasing public health burden with high mortality, highlighting the need for improved early detection strategies. Current surveillance tools, including Alpha-fetoprotein (AFP) and ultrasound, lack sufficient sensitivity for early-stage HCC detection. We analyzed serum samples from 131 patients, including 58 with cirrhosis and 73 with MASLD-related HCC (42 early-stage, 31 late-stage), using an nLC-stepped HCD-PRM-MS/MS workflow for targeted N-glycome profiling of glycopeptides derived from haptoglobin and vitronectin. Combining targeted glycopeptides with AFP significantly improved HCC detection compared with AFP alone. The optimal panel for all HCC versus cirrhosis (AFP + VTNC_169_A2G2F0S1 + VTNC_242_A3G3F2S2) achieved an AUC of 0.859 and 76.7% sensitivity at 90% specificity. For early-stage HCC, AFP + HP_184_A3G3F1S3 + VTNC_169_A2G2F0S1 yielded an AUC of 0.890 with 66.7% sensitivity at 1% specificity. A SHAP-selected Gaussian Naive Bayes model based on seven molecular/glycopeptide features, without demographic variables, further improved performance, achieving ROC-AUC values of 0.9985 in training and 1.0000 in independent testing cohorts, with accuracies of 98.1% and 100.0%, respectively.

Pediatric low-grade glioma (pLGG) is the most common type of brain tumor in children, accounting for approximately 30% of all central nervous system tumors in children. pLGG has multiple molecular subtypes that differ in disease progression, recurrence patterns, and treatment responses. Conventional wet lab approaches including molecular profiling and histopathological studies for pLGG characterization are time consuming, costly, and laborious. Recently, methods based on artificial intelligence (AI) or machine learning (ML) have been widely used for pLGG molecular categorization, but most of them can only identify two or three pLGG subtypes. To more comprehensively characterize the molecular subtypes of pLGG and their potential biological and therapeutic significance, we develop an integrated meta-clustering approach, namely Meta-pLGG, that can explore high resolution molecular subtypes and their transcriptional heterogeneity for pLGG. Specifically, we first performed multiple rounds of random projection (RP) to generate dimension-reduced feature vectors from pLGG transcriptomics data, each of which was subsequently clustered by different clustering algorithms including hierarchical clustering, K-means, Self-Organizing Maps (SOM), Non-negative Matrix Factorization (NMF), Gaussian Mixture Model (GMM), and Spectral Clustering, as base clustering methods. Then, to yield robust clustering performance, we integrated the clustering results of these RP based individual clustering algorithms by adopting a weighted meta-clustering (wMetaC) approach. Results based on 532 pLGG patients suggested that our proposed approach demonstrated superior stability and discriminative powers for higher resolution pLGG subtyping compared to conventional approaches. Based on consensus matrix analysis, we identified two major pLGG mega-subtypes, with one further subdivided into three subgroups and the other into two. Then, we performed cluster specific differential gene expression analysis, molecular pathway analysis, and gene-drug-disease association analysis. The results showed that the identified five subgroups exhibited significant subtype-specific transcriptomic heterogeneity. In summary, our meta-clustering approach demonstrated much higher performance and robustness in identifying higher resolution molecular subtypes of pLGG, revealing the molecular heterogeneity within pLGG and potentially providing new insights for more precise molecular subtyping and precision therapy.

Ancestry-related genomic variation may contribute to phenotypic diversity. However, its cellular consequences remain difficult to evaluate. To address this, we developed a human induced pluripotent stem cell (hiPSC)-based framework for investigating such phenotypes. We used the Japanese population as a model of Jomon-related and continental East Asian ancestry admixture. Jomon ancestry, estimated using the Jomon allele score (JAS), was associated with clinical measures, including renal-related traits and serum triglycerides in 34,152 Tohoku Medical Megabank participants. Peripheral blood mononuclear cell transcriptomics of donors from the upper and lower extremes of the JAS distribution (n = 10) suggested group-level differences in immune-related gene expression. To test whether the serum triglyceride association coincided with hepatic cellular phenotypes, we differentiated hiPSCs from these donors into hepatocyte-like cells. Neither intracellular triglycerides nor lipid-related gene expression differed between groups. These findings highlight the value of a genetically characterized hiPSC-based framework for investigating ancestry-associated phenotypes across cellular contexts.

We applied the latest RIC-seq technology to study the RNA-RNA interactions in human liver HepG2 cells. Integrating ChIP-seq data with the RIC-seq data, we clas-sified vari ous transcription factors (TFs) into two distinct clusters based on the in-teraction enrichment scores of the uaRNA (RNA transcribed in promoter regions) and eRNA (RNA tran scribed in enhancer regions) potentially transcribed by the TF binding sites. We went on to compare the genomic features and clinical significance across different clusters, by further integrating Hi-C, ATAC-seq, and clinical data from liver cancer patients. Our results showed that Cluster 1 exhibited low spatial interaction (RIC-seq derived RNA-RNA interactions between regulatory RNA ele-ments) and high tissue specificity (high proportion of binding sites covered by shared peaks between two different cell lines for each histone modification marker), showed higher associations with the B compartments and nucleoli, contained less densely clustered TF binding sites, and is more correlated with clinical outcomes. In contrast, Cluster 2 has higher spatial interaction, is more likely to be found in the A compartments and nuclear speckles, contained more densely clustered TF binding sites, and is more prone to form loops. These findings suggested that these transcrip-tion factors might play critical roles in regulating gene expression and serve as po-tential biomarkers and therapeutic targets with significant implications for cancer prognosis and treatment. In summary, groups of TFs could occur as distinct sub-networks based on spatial proximity of non-coding RNAs transcribed by the TF binding sites. These sub-networks appear to be related to regulatory differences for the TF groups and may partially explain the chromosome compartments/domains.

A protein whose removal or deletion causes significant disruption or collapse of a protein-protein interaction (PPI) network is referred to as a vulnerable protein. Such proteins may serve as valuable therapeutic or diagnostic targets in disease-associated networks. In this study, two PPI networks were constructed, one for HPV-positive and the other for HPV-negative head and neck squamous cell carcinoma (HNSCC), and the vulnerable proteins of these networks were identified by the node deletion approach. After analyzing the networks, 27 unique vulnerable proteins in HPV-positive and 72 unique vulnerable proteins in HPV-negative HNSCC were identified. Among them, one HPV-positive and seven HPV-negative HNSCC vulnerable proteins were further chosen by integrating multi-omics data. To exploit the vulnerabilities of these proteins, candidate synthetic lethal (SL) partners were predicted whose inhibition may selectively impair tumor survival. Subsequently, drug-gene interaction analysis was performed to identify inhibitors targeting the SL partners of these vulnerable proteins. Notably, in HPV-positive HNSCC, TOP2A, CHEK1, and CHEK2 genes were identified as SL partners of TTN, and their inhibitors were already clinically approved. While in HPV-negative HNSCC, ADA and MMP19 were identified as an SL partner of LMO7; TMEM45B, CDH3, and ELF3 genes were identified as an SL partner of CGN; and ZNF433 was identified as an SL partner of FLNC. However, MMP19, ZNF433, and TMEM45B inhibitors were not reported. Thus, these vulnerable proteins, including their SL partners, provide novel avenues to explore and develop more efficient and precise therapeutic and diagnostic strategies.

Desertification exerts dramatic selection pressures on the evolution of plants. Despite the key role of ecological adaptation by natural selection to arid grasslands and subsequent intraspecific divergence, specific mechanisms driving this process remain poorly understood. Psammochloa villosa, a perennial forage grass endemic to the arid grasslands in Northwest China, where it thrives in shifting and semi-fixed sand land due to its exceptional drought tolerance, provides an ideal system to study adaptive evolution to aridity. In our study, we assembled a high-quality, chromosome-scale genome and conducted genomic resequencing of 42 populations across its major distribution. The genome assembly, which is approximately 1.55 Gb in size, has a super-scaffold N50 of 66.79 Mb, with 75.84% of the sequences identified as transposable elements. Coalescent phylogeny and genomic collinearity analyses strongly supported that P. villosa and Neotrinia splendens, as the closest taxa, shared a recent whole-genome duplication (WGD) event occurring approximately 18-20 Mya and followed by their divergence around ~11.2 Mya. Based on ancestral grass karyotype (AGK) reconstruction from synteny analysis, our results suggest that, relative to the AGK after the {rho}-WGD event, P. villosa and N. splendens underwent similar chromosomal restructuring and lineage-specific retention of numerous copies, providing a genomic basis for potential ecological adaptation and intraspecific diversification. The expanded XTH family, which encodes enzymes mediating xyloglucan endotransglucosylation and hydrolysis and thereby regulating xyloglucan remodeling, showed strong transcriptional responses under PEG-6000 treatment, suggesting that retained copies may be associated with xerophytic adaptation in P. villosa. Together, these findings suggest that WGD-derived gene retention created a delayed reservoir of genetic diversity that was later shaped by desertification and Qinghai-Xizang Plateau environmental changes, contributing to climate-associated genomic islands and intraspecific differentiation.

Archaic introgression has shaped the evolutionary history of Eurasian populations, yet Central Eurasian region remains understudied despite being at the crossroads of ancient human migration. Here, we analyzed the whole-genome data of five Central Eurasian (CE) individuals from Early Bronze Age (EBA) and five present-day CE individuals to characterize the archaic introgression landscape. We estimated that archaic introgression from Neanderthal and Denisovan archaic hominins comprises approximately 2.2% of the Central Eurasian genomes. Both amount and chromosomal distribution of archaic introgression remained largely unchanged between the EBA and present-day CE individuals. Putative introgressed fragments matching the Altai Neanderthal and the Altai Denisovan were retrieved. Our results suggest that while the archaic introgression levels seemingly remained stable over the past several thousand years, larger modern CE genomes panels will be required to fully characterize the genomic landscape of archaic ancestry in the region.

Background: Alpha-mannosidosis (AM) is a debilitating lysosomal storage disorder (LSD) caused by insufficient activity of alpha-mannosidase, an enzyme responsible for cleavage of alpha-linked mannose residues from N-glycan-type oligosaccharides during their degradation. As a result, specific oligosaccharides accumulate in tissues and are also detectable in body fluids. AM is a progressive multi-organ disease for which hematopoietic stem cell transplantation (HSCT) and enzyme replacement therapy (ERT) are therapeutic options, however, treatment results depend considerably on early diagnosis. Therefore, inclusion of AM into newborn screening panels would significantly contribute to improving therapeutic outcome. Methods: A rapid and fully quantitative method for the specific AM oligosaccharide biomarker Man2GlcNAc in dried blood spots (DBS) using ultra-performance liquid chromatography coupled to tandem-mass spectrometry was developed. DBS from 420 controls, 27 AM patients (including five newborn screening samples), and 35 patients with other LSDs were analyzed. Results: A distinct separation in biomarker concentrations was observed between unaffected individuals and untreated AM patients (0.1 - 0.8 vs 3.2 - 11.3 mol Man2GlcNAc/L, p<0.0001). Furthermore, Man2GlcNAc concentrations significantly differed between untreated patients and those on ERT (n = 17, p<0.0001) or after HSCT (n = 15, p<0.0001), and between both treatment groups (p<0.0001). Moreover, we demonstrated that this method is appropriate for monitoring treatment efficacy by analysis of samples from individual patients before and after initiation of therapy. We retrospectively analyzed residual samples of dried blood spots (DBS) from routine newborn screening that had been stored from patients with later-confirmed AM, which resulted in clear separation from the concentrations observed in a newborn screening control cohort (7.7 - 11.3 vs 0.1 - 0.8 mol Man2GlcNAc/L, p=0.0001). Conclusion: We present a simple and rapid method for quantitative analysis of the specific AM oligosaccharide biomarker Man2GlcNAc, which is remarkably stable compared to enzymatic assays. Thus, the method is useful for selective screening in suspected cases as well as for monitoring the response to different forms of treatment. Furthermore, Man2GlcNAc has the potential to serve as a useful biomarker for newborn screening offering the prospect of earlier treatment initiation, improved therapeutic outcome and higher quality of life for AM patients.

The long-standing challenge of seasonal influenza vaccines providing poor protection is attributed to the rapid and continuous evolution of the virus. In theory, the recommendation of vaccine strain is to predict the dominant strain with best fitness in the upcoming season. Here, we develop FutureFlu, a biologically-grounded framework that quantifies the fitness score of influenza variants in upcoming seasons, by integrating metrics on three levels: molecular genetic divergence, individual immune escape, and population-scale transmission dynamics. The fitness scores of influenza variants present strong positive correlation with their actual observed frequencies in the next seasons. Validation across 24 seasons of three influenza subtypes shows FutureFlu recommends antigenically matched vaccine strains more frequently than annual recommendations, particularly for challenging subtypes: 83.3% versus 45.8% seasons for A/H3N2, and 75.0% versus 33.3% seasons for B/Victoria. Furthermore, FutureFlu significantly outperforms the currently used methods in vaccine strain selection whether or not there is available antigenic data from hemagglutination inhibition (HI) assay, which provides a valuable supplement for WHO vaccine recommendation. To support global public health implementation, an open online platform (futureflu.com.cn) has been established to offer real-time viral fitness predictions and vaccine recommendations.

Despite significant advancements in microsurgical techniques in recent years, the treatment and prognosis of craniopharyngiomas remain unsatisfactory. As a central nervous system tumor located adjacent to important brain structures such as the hypothalamus-pituitary axis and accompanied by a highly inflammatory microenvironment, the tumor heterogeneity and tumor microenvironment characteristics of papillary craniopharyngiomas (PCPs) remain unclear. In this study, we integrated multimodal single-cell and spatial profiling from PCP tissue and peripheral blood mononuclear cells (PBMCs) to elucidate the tumor heterogeneity and microenvironment characteristics of PCP. Our single-cell and spatial analyses defined four specific tumor cell states in PCP, representing specific transcriptional regulatory programs and spatial heterogeneity characteristics during tumor progression. By constructing a spatial niche composed of tumor, immune, and stromal cells, we analyzed the cellular and spatial ecosystem of PCP at multiple levels to further assess the communication relationships between different tumor cell states and microenvironment cells. This study established a multidimensional molecular atlas of PCP from the perspectives of cell state, spatial structure, and microenvironment interactions, providing a foundation for understanding its biological behavior and exploring new intervention strategies.

Motivation: MicroRNA-target interaction prediction remains challenging because many existing tools provide prediction scores or ranked candidate lists without making the supporting evidence easy to interpret or relate to a specific biological context. Results: Here, we developed miRAssist, a context-aware evidence-integration framework for interpretable miRNA-target prioritization. miRAssist integrates six evidence families, including sequence complementarity, thermodynamic stability, sequence conservation, target-site accessibility, functional binding, and functional repression. A sequence-defined candidate universe was generated, resulting in 280,917 candidate interactions. Using miRTarBase-supported interactions as known-positive labels, six supervised scoring approaches were evaluated using a grouped train/test split by miRNA. Random forest showed the strongest performance and was selected. miRAssist also produced stronger known-positive enrichment than established miRNA-target prediction models in the evaluated benchmark. An LLM-assisted interface further supports natural-language database querying and evidence-grounded summarization of prioritized candidates.

Background: Large-scale biobanks have advanced genomic and epidemiologic research, but many rely on infrequent biological sampling and limited digital phenotyping. The Yamanashi Multi-omics Cohort (YMoC) was established to support longitudinal assessment of molecular, clinical, and behavioural changes in a screening-defined cohort of adults at elevated metabolic risk without diagnosed diabetes. Methods: YMoC is a longitudinal cohort of 215 adults aged 30-70 years in Yamanashi Prefecture, Japan, who met prespecified glycaemic eligibility criteria at health check-up, including fasting plasma glucose 100-125 mg/dL (5.6-6.9 mmol/L) and HbA1c <6.5%. Participants underwent three in-person visits over six months. Measurements include 75-g oral glucose tolerance testing with serial sampling, clinical biochemistry, anthropometry, liver elastography, and collection of blood, urine, stool, and saliva for multi-omics profiling. Between visits, participants wore a Fitbit Inspire 3 and completed daily app-based questionnaires using the Taohealth app. Current molecular data include genome-wide single nucleotide polymorphism array genotyping and longitudinal plasma proteomics in a subset. Conclusions: YMoC is designed to evaluate within-person molecular and phenotypic trajectories in a screening-defined metabolic-risk cohort. The cohort provides a dense longitudinal resource linking clinical assessments, biospecimens, omics assays, and digital phenotyping, including analyses of insulin-resistance-related markers such as homeostasis model assessment of insulin resistance (HOMA-IR).

Recursive splicing (RS) is a process in which an intron is removed from a nascent RNA molecule in two or more splicing events rather than one. We introduce a novel approach for detecting recursive splice sites (RSSs), the intronic loci at which RS events occur, based on alignment of total RNA-seq data to short, customized "target" sequences. We applied this approach to a data set from a recent study of gene expression in the human brain, using parameters corresponding to a very low false discovery rate, and found 3,022 RSSs that appear in 2,775 distinct introns from 2,407 genes. 2,891 (96%) of these RSSs are in protein-coding genes. The median length of recursively spliced introns from this set is 10,114 base pairs, which is substantially longer than the median human intron, but much shorter than average RS intron lengths reported in prior studies. Our work dramatically increases the number of known RSSs in the human genome and provides a generalizable bioinformatics pipeline for annotating RSSs from total RNA-seq data.

Detecting somatic mutations is a cornerstone of cancer genomics and clinical genotyping, traditionally relying on whole exome (WES) or whole genome sequencing (WGS) data. While RNA sequencing (RNA-seq) is primarily used to assess gene expression, it can also be utilized for somatic variant detection. RNA-seq inherently captures expressed variants, thereby reducing the identification of passenger mutations and mitigating the annotation biases associated with WES. We performed somatic variant calling and cancer driver detection from tumor RNA-seq data from 1,359 RNA-seq samples with matched WGS from the Pan-Cancer Analysis of Whole Genomes (PCAWG) Consortium of the International Cancer Genome Consortium (ICGC) and The Cancer Genome Atlas (TCGA). Our approach using RNA-seq data, alone, identified at least one putative driver gene in 15 of 16 cancer types, consistent with driver predictions from WGS. We further observed that the variant allele frequency (VAF) of RNA-seq-derived variants reflects their predicted functional impact, providing information specific to transcriptomic data. Because RNA-seq data can be generated and analyzed at lower computational cost and time than WGS or WES, it represents a resource-efficient alternative for cancer driver detection. This study provides an extensive pan-cancer resource of RNA-seq-derived variants from the PCAWG consortium and supports the utility of RNA-seq as a complementary approach for characterizing expressed somatic variation in cancer.

Biomarker discovery from high-dimensional RNA sequencing data remains challenging. Conventional methods such as differential expression analysis, pairwise protein-protein interaction networks, and weighted gene co-expression network analysis suffer from false-positive interactions, transitivity-driven noise, and arbitrary thresholds for feature selection. Addressing these limitations, this study proposed a novel framework, Hypergraph Regularised Fuzzy Rough Gene Selection (HR-FRGS), for knowledge-driven biomarker discovery. This pipeline integrates information-rich protein clustering with entropy filtering and PPI-pruned co-expression interactions to form a dual-layer hypergraph. A hypergraph-based scoring method that combines Random Walk with Restart diffusion and Hypergraph Betweenness Centrality is used to identify genes that are central in the biological network. Furthermore an autonomous fuzzy rough set selection eliminates arbitrary threshold dependency for gene selection. HR-FRGS was applied to four TCGA cancer cohorts: Lung Adenocarcinoma, Head and Neck Cancer, Kidney Clear Cell Carcinoma, and Colon Cancer, which reduced over 60,000 transcripts to compact biomarker panels. Validation using classical Machine Learning algorithms with cross-validation demonstrated that biomarker panels matched or exceeded the classification performance of the full transcriptome (AUC > 0.99, MCC > 0.94) and confirmed generalizability (MCC $\geq$ 0.893). Benchmarking against baseline methods and an ablation study showed the contribution of major parts of the algorithm. Functional enrichment analysis captures established pan-cancer hallmarks and cohort-specific oncogenic mechanisms, confirming the biological relevance of the identified genes.

Short-read sequencing-based single-cell transcriptomics represents the current gold standard for studying cellular transcriptomes but remains limited in its ability to resolve full-length transcript isoforms and splicing patterns. Long-read single-cell and single-nucleus RNA sequencing (LR sc/snRNA-seq) enables the transcriptome-wide characterization of full-length isoforms at cellular resolution, yet the relative performance of commercially available workflows remains insufficiently explored. Here, using nuclei extracted from a standardized multi-species benchmark sample and Oxford Nanopore Technologies long-read sequencing, we systematically benchmarked four LR snRNA-seq strategies: 10x Genomics 3', 10x Genomics 5', ArgenTag, and Parse Biosciences. Comparing transcriptome features qualitatively and quantitatively, as well as the concordance with matched short-read data and the ability to resolve cellular heterogeneity, we identified substantial method-specific differences in read length and yield, transcript coverage, isoform detection, and recovery of sample-specific biological information, with the 10x Genomics 3' and 5' assays emerging as the most balanced approaches for comprehensive isoform-resolved single-nucleus transcriptomics. Altogether, our study provides a systematic assessment of four commercially available workflows for performing LR snRNA-seq and highlights key methodological trade-offs related to distinct library preparation strategies, thus providing practical guidance for future isoform-resolved transcriptome studies at the single-nucleus level.

Glioblastoma (GBM) is a highly aggressive brain tumor with an extremely poor 5-year survival rate of 6.9%, largely attributable to the lack of reliable biomarkers. While competing endogenous RNA (ceRNA) and copy number variation (CNV) analyses offer unique biomarker identification potential, current approaches neglect the integration of multiple regulatory mechanisms for biomarker detection. To address this limitation, we applied relational graph convolutional networks (RGCNs) to ceRNA and CNV knowledge graphs through a novel late fusion ensemble architecture. The proposed architecture outperformed baseline models and identified five novel biomarkers, including hsa-miR-196a and hsa-miR-224. Kaplan-Meier survival analysis and Cox regression indicated that the identified genes hold significant prognostic and diagnostic power. The early stratification of the Kaplan-Meier curves indicates the potential these genes hold for patient survival prediction. The results illustrate that a late fusion RGCN ensemble effectively captures complex gene interactions, overcoming limitations of existing models and providing a framework for biomarker discovery. The novel biomarkers serve as prospective targets for future GBM therapeutic development and candidates for non-invasive diagnostic assays.

The attention mechanism offers the possibility for data-driven discovery of biological principles. However, for important protein families such as human olfactory receptors, the extent to which attention can associate with biologically meaningful key regions lacks systematic validation. In this study, using human olfactory receptors (ORs) as a model, we constructed CrossVOI, a VOC-OR interaction prediction framework based on protein language models and cross-attention, achieving predictive performance superior to existing methods. Furthermore, we systematically analyzed the attention distributions of CrossVOI and found that attention not only focused on ligand-binding interfaces and evolutionarily conserved sites, but also to some extent identified certain dynamically regulated regions. In summary, we propose CrossVOI, currently the best-performing framework for VOC-OR interaction prediction, and analyze the interpretability of the attention mechanism for human ORs. This study provides insights into the interpretability of protein function prediction methods and is expected to contribute to the exploration of attention mechanisms in biological mechanisms, and provide assistance for large-scale screening and mechanistic analysis of olfactory receptors.

N6-methyladenosine (m6A) modification is the most predominant and ubiquitous internal modification of RNA in eukaryotes, serving as a key post-transcriptional regulator of gene expression that is dynamically modulated by methyltransferases (writers) and demethylases (erasers). However, while the functions of m6A methylases have been partially elucidated in insects, the identity of m6A erasers in arthropods and their chemical catalytic mechanisms, as well as biological functions, remains largely enigmatic. Here, we uncovered 2499 putative methylase genes and 1148 putative demethylase genes in 266 insect genomes, and demonstrated that ALKBH4 functions as an m6A demethylase in the whitefly, Bemisia tabaci, catalyzing the oxidative reversal of mRNA m6A modifications both in vitro and in vivo. Furthermore, we established that ALKBH4, in coordination with other core components of the m6A pathway, fulfills an essential function in regulating the transcript stability of Imaginal Disk Growth Factor 1 (IDGF1) during whitefly development. Collectively, our findings expand the evolutionary scope of the eukaryotic m6A modification system, and reveal a conserved yet insect-specific epitranscriptomic regulatory mechanism governing fundamental physiological processes and adaptive phenotypes. Significance statementThe addition of a methyl group to the N6-position of adenosine (m6A) is a highly abundant chemical modification of RNA. However, the functional role of m6A in insects and the key enzymes that regulate its levels remains poorly understood. In this study, we explored putative methylase genes and demethylase genes in hundreds insect genomes, and identified an m6A RNA demethylase, ALKBH4, in the whitefly, Bemisia tabaci. We demonstrate that ALKBH4 oxidatively reverses mRNA methylation in vivo and in vitro, in combination with other components of the m6A pathway, plays an important role in whitefly development. These findings provide new insight into m6A methylation system of insect.